﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace SmartMathLibrary.NumericalIntegration
{
    /// <summary>
    /// This class provides bivariate numerical integration by using a dynamic grid. EXPERIMANTEL CODE. 
    /// </summary>
    [Serializable]
    public class BivariateDynamicGridIntegrator : AbstractMultivariateIntegrator
    {
        /// <summary>
        /// This field holds the first dimension lower bound of the integral.
        /// </summary>
        private double a1;

        /// <summary>
        /// This field holds the second dimension lower bound of the integral.
        /// </summary>
        private double a2;

        /// <summary>
        /// This field holds the first dimension upper bound of the integral.
        /// </summary>
        private double b1;

        /// <summary>
        /// This field holds the second dimension upper bound of the integral.
        /// </summary>
        private double b2;

        /// <summary>
        /// This field holds the bivariate integral.
        /// </summary>
        private BivariateIntegral integral;

        /// <summary>
        /// This field holds the grid integrals.
        /// </summary>
        private BivariateIntegral[,] splittedIntegrals;

        /// <summary>
        /// Initializes a new instance of the <see cref="BivariateDynamicGridIntegrator"/> class.
        /// </summary>
        /// <param name="integral">The bivariate integral to integrate.</param>
        public BivariateDynamicGridIntegrator(BivariateIntegral integral)
        {
            if (integral == (BivariateIntegral) null)
            {
                throw new ArgumentNullException("integral");
            }

            this.integral = integral;
        }

        /// <summary>
        /// Gets or sets the bivariate integral to integrate.
        /// </summary>
        /// <value>The bivariate integral to integrate.</value>
        public BivariateIntegral Integral
        {
            get { return integral; }
            set { integral = value; }
        }

        /// <summary>
        /// Integrates the integral in a numerical way. After running the numerical integration
        /// it should be checked if the result has a precision error.
        /// </summary>
        /// <param name="iterations">The maximum numbers of iterations.</param>
        /// <param name="precision">The minimum precision of the numerical integration.</param>
        /// <returns>The area of the specified integral.</returns>
        public override double Integrate(int iterations, double precision)
        {
            return this.Integrate(10, 10, 20, 20, iterations, precision);
        }

        /// <summary>
        /// Integrates the integral in a numerical way. After running the numerical integration
        /// it should be checked if the result has a precision error.
        /// </summary>
        /// <param name="gridX">The number of integrals in x direction.</param>
        /// <param name="gridY">The number of integrals in y direction.</param>
        /// <param name="gridPartsX">The number of regtangles for each grid item in x direction (start value).</param>
        /// <param name="gridPartsY">The number of regtangles for each grid item in y direction (start value).</param>
        /// <param name="iterations">The maximum numbers of iterations.</param>
        /// <param name="precision">The minimum precision of the numerical integration.</param>
        /// <returns>The area of the specified integral.</returns>
        public double Integrate(int gridX, int gridY, int gridPartsX, int gridPartsY, int iterations, double precision)
        {
            if (gridX == 1)
            {
                gridX++;
            }

            if (gridY == 1)
            {
                gridY++;
            }

            double sign1 = 1;
            double sign2 = 1;
            this.a1 = this.integral.A1;
            this.a2 = this.integral.A2;
            this.b1 = this.integral.B1;
            this.b2 = this.integral.B2;

            if (b1 < a1)
            {
                double tempuri = b1;

                b1 = a1;
                a1 = tempuri;
                sign1 = -1;
            }

            if (b2 < a2)
            {
                double tempuri = b2;

                b2 = a2;
                a2 = tempuri;
                sign2 = -1;
            }

            int i = 0;

            if (ExtendedMath.IsOdd(gridPartsX))
            {
                gridPartsX++;
            }

            if (ExtendedMath.IsOdd(gridPartsY))
            {
                gridPartsY++;
            }

            gridPartsX /= 2;
            gridPartsY /= 2;

            this.SplitIntegral(gridX, gridY);
            gridX--;
            gridY--;

            int[,] gridPartMapX = new int[gridX,gridY];
            int[,] gridPartMapY = new int[gridX,gridY];
            double[,] oldGridResults = new double[gridX,gridY];
            double[,] newGridResults = new double[gridX,gridY];
            InternalBivariateGridIntegrator[,] gridIntegrators = new InternalBivariateGridIntegrator[gridX,gridY];

            for (int j = 0; j < gridX; j++)
            {
                for (int l = 0; l < gridY; l++)
                {
                    gridPartMapX[j, l] = gridPartsX;
                    gridPartMapY[j, l] = gridPartsY;
                    gridIntegrators[j, l] = new InternalBivariateGridIntegrator(this.splittedIntegrals[j, l]);
                    oldGridResults[j, l] = gridIntegrators[j, l].Integrate(gridPartsX, gridPartsY);
                }
            }

            bool[,] precisionMap = new bool[gridX,gridY];

            while (i < iterations)
            {
                for (int j = 0; j < gridX; j++)
                {
                    for (int l = 0; l < gridY; l++)
                    {
                        if (!precisionMap[j, l])
                        {
                            gridPartMapX[j, l] *= 2;
                            gridPartMapY[j, l] *= 2;
                            newGridResults[j, l] = gridIntegrators[j, l].Integrate(gridPartMapX[j, l],
                                                                                   gridPartMapY[j, l]);
                        }
                        else
                        {
                            newGridResults[j, l] = oldGridResults[j, l];
                        }
                    }
                }

                precisionMap = CreatePrecisionMap(oldGridResults, newGridResults, precision);

                if (CheckPrecision(oldGridResults, newGridResults, precision))
                {
                    this.neededIterations = i;
                    this.precisionError = false;
                    this.relativeError = GetBiggestError(oldGridResults, newGridResults);

                    if (sign1 == sign2)
                    {
                        return ComputeSum(newGridResults);
                    }

                    return -ComputeSum(newGridResults);
                }

                oldGridResults = (double[,]) newGridResults.Clone();
                i++;
            }

            this.relativeError = GetBiggestError(oldGridResults, newGridResults);
            this.neededIterations = i;
            this.precisionError = true;

            if (sign1 == sign2)
            {
                return ComputeSum(newGridResults);
            }

            return -ComputeSum(newGridResults);
        }

        private void SplitIntegral(int firstDimension, int secondDimension)
        {
            double xStart = Math.Min(a1, b1);
            double yStart = Math.Min(a2, b2);
            double xDirection = (Math.Max(a1, b1) - Math.Min(a1, b1)) / (firstDimension - 1);
            double yDirection = (Math.Max(a2, b2) - Math.Min(a2, b2)) / (secondDimension - 1);

            this.splittedIntegrals = new BivariateIntegral[firstDimension - 1,secondDimension - 1];

            for (int x = 0; x < firstDimension - 1; x++)
            {
                for (int y = 0; y < secondDimension - 1; y++)
                {
                    this.splittedIntegrals[x, y] = new BivariateIntegral(xStart + (x * xDirection),
                                                                         yStart + (y * yDirection),
                                                                         xStart + (x + 1) * xDirection,
                                                                         yStart + (y + 1) * yDirection,
                                                                         this.integral.Function);
                }
            }
        }

        private static bool CheckPrecision(double[,] a, double[,] b, double precision)
        {
            for (int i = 0; i < a.GetLength(0); i++)
            {
                for (int j = 0; j < a.GetLength(1); j++)
                {
                    if (Math.Abs(a[i, j] - b[i, j]) > precision)
                    {
                        return false;
                    }
                }
            }

            return true;
        }

        private static double GetBiggestError(double[,] a, double[,] b)
        {
            double result = double.MinValue;

            for (int i = 0; i < a.GetLength(0); i++)
            {
                for (int j = 0; j < a.GetLength(1); j++)
                {
                    if (Math.Abs(a[i, j] - b[i, j]) > result)
                    {
                        result = Math.Abs(a[i, j] - b[i, j]);
                    }
                }
            }

            return result;
        }

        private static bool[,] CreatePrecisionMap(double[,] a, double[,] b, double precision)
        {
            bool[,] result = new bool[a.GetLength(0),a.GetLength(1)];

            for (int i = 0; i < a.GetLength(0); i++)
            {
                for (int j = 0; j < a.GetLength(1); j++)
                {
                    if (Math.Abs(a[i, j] - b[i, j]) > precision)
                    {
                        result[i, j] = false;
                    }
                    else
                    {
                        result[i, j] = true;
                    }
                }
            }

            return result;
        }

        private static double GetBiggestNumber(double[,] a)
        {
            double result = double.MinValue;

            for (int i = 0; i < a.GetLength(0); i++)
            {
                for (int j = 0; j < a.GetLength(1); j++)
                {
                    if (Math.Abs(a[i, j]) > result)
                    {
                        result = Math.Abs(a[i, j]);
                    }
                }
            }

            return result;
        }

        private static double ComputeSum(double[,] a)
        {
            double result = 0;

            for (int i = 0; i < a.GetLength(0); i++)
            {
                for (int j = 0; j < a.GetLength(1); j++)
                {
                    result += a[i, j];
                }
            }

            return result;
        }
    }
}